#include "manager.h"

#ifdef _HAVE_PETSC_

int unifyVolumesAmongProcessors(pMesh theMesh, const set<int> &setOfDomains, const string &whatUnify)
{
	if (P_size() == 1) return 0;

	PetscErrorCode ierr;
	int i, col = 0;

	Mat exchangeVolumes; // matrix to store nodes from all processes

	set<int>::iterator SIter;
	for (SIter = setOfDomains.begin(); SIter != setOfDomains.end(); SIter++)
	{
		const int dom = *SIter;	// domain flag

		// parallel matrix to unify volumes from nodes with remote copies
		MatCreateMPIAIJ(PETSC_COMM_WORLD, PETSC_DECIDE,PETSC_DECIDE, M_getMaxVId(theMesh),
				1,0,PETSC_NULL,0,PETSC_NULL,&exchangeVolumes);

		pEntity node;
		list<int> remoteNodeIds;

		// loop over nodes to get those with remote copies
		VIter vit = M_vertexIter(theMesh);
		while (node = VIter_next(vit))
		{
			if (M_numRemoteCopies(theMesh,node))
			{
				i = EN_id(node)-1;
				remoteNodeIds.push_back(i);
				double val = .0;
				char volString[256];

				if ( !whatUnify.compare("volume") )
					sprintf(volString,"%s-%d",whatUnify.c_str(),dom);
				else
					sprintf(volString,"%s",whatUnify.c_str());				

				EN_getDataDbl(node,MD_lookupMeshDataId(volString),&val);
				ierr = MatSetValue(exchangeVolumes,i,col,val,ADD_VALUES); CHKERRQ(ierr);
			}
		}
		VIter_delete(vit);
		ierr = MatAssemblyBegin(exchangeVolumes,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
		ierr = MatAssemblyEnd(exchangeVolumes,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
		//ierr = MatView(exchangeVolumes,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

		// transfer all remote nodes from a set container to an array
		int numNodes = remoteNodeIds.size();
		int *nodeIDs = new int[numNodes];

		list<int>::iterator Iter;
		for (i = 0, Iter = remoteNodeIds.begin(); Iter != remoteNodeIds.end(); Iter++, i++) nodeIDs[i] = *Iter;
		remoteNodeIds.clear();

		// indices used for petsc function    
		IS isrow, iscol;
		ISCreateGeneral(PETSC_COMM_WORLD,numNodes,nodeIDs,&isrow);
		ISCreateGeneral(PETSC_COMM_WORLD,1,&col,&iscol);
		delete[] nodeIDs; nodeIDs = 0;

		// now, all nodes, regardless their location(processor), are updated.
		Mat unifiedVolumes; // matrix with all requested nodes
		ierr = MatGetSubMatrix(exchangeVolumes, isrow, iscol, PETSC_DECIDE, MAT_INITIAL_MATRIX,	&unifiedVolumes);CHKERRQ(ierr);

		ierr = MatDestroy(exchangeVolumes); CHKERRQ(ierr);
		ISDestroy(isrow);
		ISDestroy(iscol);

		// each process must get their volumes from the remote nodes
		// now, the control volume from all nodes have the same values as they were calculated in a sequential
		// manner
		PetscInt m, n;
		ierr = MatGetOwnershipRange(unifiedVolumes,&m,&n);

		int N = n - m; // number of rows in unifiedVolumes owned to each process
		double *volumes = new double[N];
		int *rowIndices = new int[N];

		for (i=0; i<N; i++) rowIndices[i] = m + i;

		/// VOLUMES 
		// transfer unified volumes to an array
		ierr = MatGetValues(unifiedVolumes,N,rowIndices,1,&col,volumes);CHKERRQ(ierr);
		ierr = MatDestroy(unifiedVolumes); CHKERRQ(ierr);

		// update mesh data structure with the unified volumes
		i = 0;
		vit = M_vertexIter(theMesh);
		while (node = VIter_next(vit))
		{
			if (M_numRemoteCopies(theMesh,node))
			{
				char volString[256];
				if ( !whatUnify.compare("volume") )
					sprintf(volString,"%s-%d",whatUnify.c_str(),dom);
				else
					sprintf(volString,"%s",whatUnify.c_str());
				EN_attachDataDbl(node,MD_lookupMeshDataId(volString),volumes[i++]);
			}
		}
		VIter_delete(vit);

		delete[] rowIndices; rowIndices = 0;
		delete[] volumes; volumes = 0;
	}
	return 0;
}
#endif //-D_HAVE_PETSC_
